\chapter{随机差分方程的前向求解}
对于如下差分方程，
\[ y_t=\lambda y_{t+1}+\gamma x_t,\qquad |\lambda|<1, t\ge 0 \]

Sargent (1987)证明了它存在唯一解。查分方程的解指的是存在一系列的值满足差分方程，比如上式就是要找一系列的$ y_t $，可以把它简记为$ (y_t)^\infty_t=0 $，它要满足上述差分方程。

要找到这个$ (y_t)^\infty_t=0 $，可以通过前向不断迭代得到该方程的解，初值任意选，比如$ (0)^\infty_t=0 $，那么依次迭代就有，
\begin{align*}
	y_t^1=&\gamma x_t,\qquad t\ge 0\\
		y_t^2=&\gamma x_t+\gamma\lambda x_{t+1},\qquad t\ge 0\\
		y_t^2=&\gamma x_t+\gamma\lambda x_{t+1}+\gamma\lambda^2x_{t+2},\qquad t\ge 0\\
		\vdots&\\
				y_t^n=&\gamma\sum_{i=0}^{n-1}\lambda^i x_{t+i},\qquad t\ge 0
	\end{align*}



\chapter{差分方程的解的性质}\label{chp_de}
\section{差分方程解的性质一般讨论}
差分方程可以写成，
\[ \Delta\bm{X}_{t+1}=\bm{AX}_t+\bm{b} \]
其中，$ \bm{X} $是$ n\times 1 $的向量，$ \bm{A} $是$ n\times n $的矩阵，$ \bm{b} $是$ n\times 1 $的向量。展开后可以写成，
\[ \bm{X}_{t+1}=\bm{\hat{A}X}_t+\bm{b} \]
其中，$ \bm{\hat{A}}=\bm{A}+\bm{I} $。很容易看到这个方程的稳态解为，
\begin{align*}
	\bm{X}_{t+1}=&(\bm{A}+\bm{I})\bm{X}_t+\bm{b}\\
	\Longrightarrow 	\bar{\bm{X}}=&\bm{A}\bar{\bm{X}}+\bar{\bm{X}}+\bm{b}\\
	\Longrightarrow \bar{\bm{X}} =& -\bm{A}^{-1}\bm{b}
\end{align*}

继续回到原来的差分方程，若令初值为$ \bm{X}_0 $，则方程可以后向递归得到解，
\begin{align}\nonumber
	\bm{X}_1=&\bm{\hat{A}X}_0+\bm{b}=(\bm{A+I})\bm{X}_0+\bm{b}\\\nonumber
	\bm{X}_2=&\bm{\hat{A}X}_1+\bm{b}=\bm{\hat{A}^2X}_0+\bm{\hat{A}b}+\bm{b}=(\bm{A+I})^2\bm{X}_0+(\bm{A}+\bm{I})\bm{b}+\bm{b}\\\nonumber
	\vdots&\\\label{de_eq_1}
	\bm{X}_t=&\bm{\hat{A}X}_{t-1}+\bm{b}=\bm{\hat{A}}^t\bm{X}_0+\sum_{i=0}^{t-1}\bm{\hat{A}}^{t-i-1}\bm{b}=(\bm{A+I})^t\bm{X}_0+\sum_{i=0}^{t-1}(\bm{A+I})^{t-i-1}\bm{b}
\end{align}

可以看到，只要给定初值，就可以把$\bm{X}_t,t\in 1,2,\cdots  $计算出来，而且关键就是计算矩阵的幂。矩阵的幂，只要做一个特征值分解，矩阵的幂是非常容易计算的。矩阵$ \bm{A} $如果是可对角化的(即便不可以对角化，还可以使用Jordan分解)，那么它一定有一个特征值分解，即
\[ \bm{A}= \bm{P}^{-1}\bm{\Lambda P}\]
其中$ \bm{P} $是由$ \bm{A} $的特征向量$ \bm{v}_1,\bm{v}_2,\cdots,\bm{v}_n $构成的矩阵，$ \bm{\Lambda} $是由特征值$ \lambda_1,\lambda_2,\cdots,\lambda_n $构成的对角矩阵。那么，很容易可以看到，
\[ \bm{A}^k= \bm{P}^{-1}\bm{\Lambda}^k\bm{P} \]

\eqref{de_eq_1}式可以重新写成，
\[ \bm{X}_t=\bm{\hat{A}}^t\bm{X}_0+(\bm{I}+\bm{\hat{A}}+\bm{\hat{A}}^2+\cdots+\bm{\hat{A}}^{t-1})\bm{b} \]
 因为
 \[ (\bm{I}+\bm{\hat{A}}+\bm{\hat{A}}^2+\cdots+\bm{\hat{A}}^{t-1})(\bm{I}-\bm{\hat{A}})=\bm{I}-\bm{\hat{A}}^t \]
 因此，有，
 \begin{align*}
 	\bm{X}_t=&\bm{\hat{A}}^t\bm{X}_0+(\bm{I}-\bm{\hat{A}}^t)(\bm{I}-\bm{\hat{A}})^{-1}\bm{b}\\
=&\bm{\hat{A}}^t\bm{X}_0+(\bm{I}-\bm{\hat{A}}^t)(\bm{I}-\bm{I}-\bm{A})^{-1}\bm{b}\\
=&\bm{\hat{A}}^t\bm{X}_0+(\bm{I}-\bm{\hat{A}}^t)(-\bm{A}^{-1}\bm{b})\\
=&\bm{\hat{A}}^t\bm{X}_0+(\bm{I}-\bm{\hat{A}}^t)\bar{\bm{X}}\\
=&\bm{\hat{A}}^t(\bm{X}_0-\bar{\bm{X}})+\bar{\bm{X}}
 \end{align*}

 可见，只要随着时间的增长，$ \hat{\bm{A}}^t $趋近于0，$ \bm{X}_t $就会接近稳态。现在就是看$ \hat{\bm{A}}^t $趋近于0需要满足什么条件。对$ \hat{\bm{A}}^t $进行特征值分解，有，
 \begin{equation}\label{de_eq_x}
 \bm{X}_t= \bm{P}\begin{bmatrix}
 	\hat{\lambda}^t_1&\cdots&0\\
 	\vdots&\ddots&\vdots\\
 	0&\cdots&\hat{\lambda}_n^t
 \end{bmatrix}\bm{P}^{-1}(\bm{X}_0-\bar{\bm{X}})+\bar{\bm{X}}	
 \end{equation}


如果令
\[ \bm{P}^{-1}(\bm{X}_0-\bar{\bm{X}})=\begin{bmatrix}
	c_1\\\vdots\\c_n
\end{bmatrix} \]

那么\eqref{de_eq_x}式可以写成，
\begin{equation}\label{de_eq_x2}
\bm{X}_t=c_1\bm{v}_1\hat{\lambda}^t_1+c_2\bm{v}_2\hat{\lambda}^t_2+\cdots+c_n\bm{v}_n\hat{\lambda}^t_n+\bar{\bm{X}} 	
\end{equation}


其中$ \bm{v}_i,i=1,2,\cdots,n $是特征向量，即$ \bm{P} $中的列。因此，稳态就等价于特征值的模(如果是实数，就是绝对值)要小于1，即$ |\hat{\lambda}_i|<1,i=1,2,\cdots,n $。

\section{二元情况的一个深入分析}
为了详细讨论解的性质，现在看一个二元系统，
\[\begin{bmatrix}
\Delta X_{1t}\\ \Delta X_{2t}
\end{bmatrix}=\begin{bmatrix}
a_{11}&a_{12}\\
a_{21}&a_{22}
\end{bmatrix}\begin{bmatrix}
{X}_{1t}\\ {X}_{2t}
\end{bmatrix}+\begin{bmatrix}
b_{1}\\ b_{2}
\end{bmatrix},\qquad\text{其中,}\bm{A}=\begin{bmatrix}
a_{11}&a_{12}\\
a_{21}&a_{22}
\end{bmatrix}\]

那么，给定初值$ \bm{X}_0 $，方程的解就可以写成，
\[ \bm{X}_t=(\bm{A}+\bm{I})^t(\bm{X}_0-\bar{\bm{X}})+\bar{\bm{X}} \]

根据前面的分析，方程稳不稳定依赖于$ \bm{A}+\bm{I} $(或者说$ \bm{A} $)\footnote{$ \bm{A} $的特征值为$ \lambda $，则$ \bm{A}+\bm{I} $的特征值为$ \lambda+1 $，且它们的特征向量是一样的。}的特征值。而它的特征值无非就是解方程$ |\bm{A}-\lambda\bm{I}|=0 $，这个方程的根为，
\[ \lambda_{1},\lambda_{2}=\frac{a_{11}+a_{22}\pm \sqrt{(a_{11}+a_{22})^2-4(a_{11}a_{22}-a_{12}a_{21})}}{2} \]

若定义$ D= (a_{11}+a_{22})^2-4(a_{11}a_{22}-a_{12}a_{21})$，此时就可以分三种情况讨论解的性质：

\paragraph{(1) D>0：}此时方程必然有实根，而且是不同的实根。因此，矩阵$ (\bm{A}+\bm{I}) $就是可以对角化的，也就是说可以写成，
\[ \begin{bmatrix}
	{X}_{1t}\\ {X}_{2t}
\end{bmatrix}=\bm{P}\begin{bmatrix}
(\lambda_{1}+1)^t&0\\ 0 & (\lambda_{2}+1)^t
\end{bmatrix}\bm{P}^{-1}\begin{bmatrix}
X_{1,0}-\bar{X}_1\\ X_{2,0}-\bar{X}_2
\end{bmatrix}+\begin{bmatrix}
\bar{X}_1\\\bar{X}_2
\end{bmatrix} \]

\begin{itemize}
	\item 如果$ \lambda_{1}+1<1,\lambda_{2}+1<1 $，则意味着系统会趋于稳态解。
	\item 如果$ \lambda_{1}+1>1,\lambda_{2}+1>1 $，则意味着系统会发散，因为随着时间的流逝，它远离了稳态$ \bar{\bm{X}} $。
	\item 如果$ |\lambda_{1}+1|<1,|\lambda_{2}+1|>1 $,则随着时间的流逝，$ (\lambda_{1}+1)^t $会趋于0，但$ (\lambda_{2}+1)^t $会变成无穷大。这意味着，只有当初值满足$X_{2,0}=\bar{X}_2,X_{1,0}\ne \bar{X}_1  $时，系统就消去了不稳定项，然后趋于稳态，否则系统永远是发散的。这就是常说的鞍点稳定。
	\item 如果$ |\lambda_{1}+1|=1,|\lambda_{2}+1|<1 $，则系统不会移动，当然也不会趋近稳态。另外，因为是绝对值等于1，所以当$ (\lambda_{2}+1)^t $趋于0时，系统可能会在1和-1间震荡。这个情形通常称为边际稳定(marginally stable)。当然，如果初始条件$ X_{1,0}= \bar{X}_1 $就能消去对应的项，系统又能趋于稳态。
\end{itemize}

\paragraph{(2) D=0:}此时方程有2个相等的实根，然后要看矩阵$ \bm{A} $是不是对角阵：

\begin{itemize}
	\item 如果矩阵$ \bm{A} $是对角阵，那么
	\[ \begin{bmatrix}
		\Delta{X}_{1t}\\ \Delta{X}_{2t}
	\end{bmatrix}=\begin{bmatrix}
		a_{11}&0\\ 0 & a_{22}
	\end{bmatrix}\begin{bmatrix}
		X_{1,t}\\ X_{2,t}
	\end{bmatrix}+\begin{bmatrix}
		b_1\\b_2
	\end{bmatrix} \]
	此时，$\lambda=a_{11}=a_{22} $，从而有，
\[ \begin{bmatrix}
	{X}_{1t}\\ {X}_{2t}
\end{bmatrix}=\begin{bmatrix}
	(\lambda_{1}+1)^t&0\\ 0 & (\lambda_{2}+1)^t
\end{bmatrix}\begin{bmatrix}
	X_{1,0}-\bar{X}_1\\ X_{2,0}-\bar{X}_2
\end{bmatrix}+\begin{bmatrix}
	\bar{X}_1\\\bar{X}_2
\end{bmatrix} \]
此时，如果$ |\lambda+1|<1 $	则稳定，如果$ |\lambda+1|>1 $	则发散，如果$ |\lambda+1|=1 $	则系统不会动。
	\item 如果$ \bm{A} $不是对角阵，那么要进行Jordan分解，
	\[ \begin{bmatrix}
		{X}_{1t}\\ {X}_{2t}
	\end{bmatrix}=\bm{Q}\begin{bmatrix}
		\lambda^t&t\cdot\lambda^{t-1}\\ 0 & \lambda^t
	\end{bmatrix}\bm{Q}^{-1}\begin{bmatrix}
		X_{1,0}-\bar{X}_1\\ X_{2,0}-\bar{X}_2
	\end{bmatrix}+\begin{bmatrix}
		\bar{X}_1\\\bar{X}_2
	\end{bmatrix} \]
	其中$ \lambda $是$ \bm{A}+\bm{I} $的特征根，$ \bm{Q} $是对应的特征向量。此时，如果$ |\lambda|<1 $	则稳定，如果$ |\lambda|>1 $	则发散，如果$ |\lambda|=1 $	则也不稳定了。因为此时，把Jordan分解的三个矩阵按照\eqref{de_eq_x2}式形式展开，会是$ c_1\bm{v}_1\lambda^t+c_1c_2\bm{v}_2t\cdot\lambda^{t-1} $，即便$ \lambda=1 $也会发散。
\end{itemize}

\paragraph{(3) D<0:} 此时方程仅有复根，$ \lambda +1= \alpha\pm\beta i $，但值不一样，所以可以对角化，
\[ \begin{bmatrix}
	{X}_{1t}\\ {X}_{2t}
\end{bmatrix}=\bm{P}\begin{bmatrix}
	(\lambda_{1}+1)^t&0\\ 0 & (\lambda_{2}+1)^t
\end{bmatrix}\bm{P}^{-1}\begin{bmatrix}
	X_{1,0}-\bar{X}_1\\ X_{2,0}-\bar{X}_2
\end{bmatrix}+\begin{bmatrix}
	\bar{X}_1\\\bar{X}_2
\end{bmatrix} \]
因为复数与三角函数的对应关系\footnote{这里回忆一下欧拉公式：指数函数与任意复数之间的关系。根据泰勒展开公式，有，
\[ e^x=\sum_{k=0}^{\infty}\frac{x^k}{k!},\;\; \cos x=\sum_{k=0}^\infty (-1)^k\frac{x^{2k}}{(2k)!},\;\; \sin x =\sum_{k=0}^\infty (-1)^k\frac{x^{2k+1}}{(2k+1)} \]
因此，有，
\begin{align*}
	e^{ix}=& \sum_{k=0}^{\infty}\frac{(ix)^k}{k!}=\sum_{k=0}^\infty (-1)^k\frac{x^{2k}}{(2k)!}+i\sum_{k=0}^\infty (-1)^k\frac{x^{2k+1}}{(2k+1)}\\
	=&\cos x+i\sin x
\end{align*}
进一步可知，
\[ re^{ix} =r\cos x+i\cdot r\sin x\]
上式右边可以表示任何一个复数，左边表示任何一个复数都可以写成指数形式。}，所以有，
\[ \begin{bmatrix}
	{X}_{1t}\\ {X}_{2t}
\end{bmatrix}=\bm{P}\begin{bmatrix}
	\rho^t(\cos \theta t+i \sin\theta t)&0\\ 0 & \rho^t(\cos \theta t-i \sin\theta t)
\end{bmatrix}\bm{P}^{-1}\begin{bmatrix}
	X_{1,0}-\bar{X}_1\\ X_{2,0}-\bar{X}_2
\end{bmatrix}+\begin{bmatrix}
	\bar{X}_1\\\bar{X}_2
\end{bmatrix} \]
其中，$ \rho=\sqrt{(\alpha+1)^2+\beta^2} $是复数$ \lambda+1 $的模，$ \theta=\arctan((\alpha+1)/\beta) $是它的角。
\begin{itemize}
	\item 如果$ \rho<1 $，则随着时间的流逝，系统会震荡而趋于稳定。
	\item 如果$ \rho>1 $，则随着时间的流逝，系统会震荡而趋于发散。
	\item 如果$ \rho=1 $，则系统边际稳定。举例而言，就是解的轨迹会循环画圈圈。
\end{itemize}

文件\lstinline|de_plot.mw|有一个军备竞赛的关于相位图的绘制代码。